clear;
WE = 4;SN = 2;lamda = 1852;
WE = WE * lamda;
SN = SN * lamda;
cover = 0.1;
theta = 2 * pi/3;
alpha = 1.5 * pi/180;
D0 = 110 - WE/2 * tan(alpha);
angle1 = pi/2 + alpha - theta/2;
angle2 = pi/2 - alpha - theta/2;
xlast = D0/(sin(angle1)/cos(alpha)-sin(theta/2)*tan(alpha));
pos = [xlast];
while 1
    xnew = (xlast * (1/cos(alpha) + sin(theta/2) * tan(alpha)/sin(angle2)) + D0 * (1-cover) * sin(theta/2) * (1/sin(angle2) + 1/sin(angle1)))/(cover * tan(alpha) * sin(theta/2) * (1/sin(angle2) + 1/sin(angle1)) - tan(alpha) * sin(theta/2)/sin(angle1) + 1/cos(alpha));
    if xnew > WE
        break;
    end
    pos = [pos, xnew];
    xlast = xnew;
end

%% 可视化
plot(1:length(pos), pos);
figure;
for i = 1:length(pos)
    dep = -(pos(i) * tan(alpha) + D0);
    plot3([pos(i),pos(i)],[0,SN],[dep,dep]);
    hold on;
end
xlabel('x/m');
ylabel('y/m');
zlabel('z/m');
hold off;
h = pos * tan(alpha) + D0;
len2 = h * sin(theta/2)/(sin(pi/2-alpha-theta/2));
len1 = h * sin(theta/2)/(sin(pi/2+alpha-theta/2));
x1 = pos - len1 * cos(alpha);
y1 = -h + len1 * sin(alpha);
x2 = pos + len2 * cos(alpha);
y2 = -h - len2 * sin(alpha);
figure
for i = 1:length(pos)
    if rem(i,4) == 0
        color = 'r';
    elseif rem(i,4) == 1
        color = 'b';
    elseif rem(i,4) == 2
        color = 'g';
    else
        color = 'k';
    end
    plot([pos(i),x1(i)],[0,y1(i)],color);
    hold on;
    plot([pos(i),x2(i)],[0,y2(i)],color);
    hold on;
    plot([x1(i),x2(i)],[y1(i),y2(i)],color);
    hold on;
end
xlabel('x/m')
ylabel('z/m')

